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ABSTRACT 

We present an extension of the harmonic-space maximum-entropy component 
separation method (MEM) for multi-frequency CMB observations that aUows one 
to perform the separation with more plausible assumptions about the receiver noise 
and foreground astrophysical components. Component separation is considered in the 
presence of spatially- varying noise variance and spectral properties of the foreground 
components. It is shown that, if not taken properly into account, the presence of 
spatially- varying foreground spectra, in particular, can severely reduce the accuracy 
of the component separation. Nevertheless, by extending the basic method to ac- 
commodate such behaviour and the presence of anisotropic noise, we find that the 
accuracy of the component separation can be improved to a level comparable with 
previous investigations in which these effects were not present. 
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1 INTRODUCTION 

An important stage in the reduction of CMB anisotropy 
data is the separation of the astrophysical and cosmological 
components. Several techniques have been suggested, includ- 
ing blind (Baccigalupi et al. 2000, Maino et al. 2002) and 
non-blind (Hobson et al. 1998, Bouchet and Gispert 1999, 
Stolyarov et al. 2002) approaches. Non-blind methods, such 
as the maximum-entropy method (MEM) or Wiener filter- 
ing, allow one to use all available prior information about 
the components in the separation process. 

A detailed description of the harmonic-space MEM ap- 
proach for flat patches of the sky was described by Hobson 
et al. (1998), and was extended later to the sphere by Stol- 
yarov et al. (2002; hereinafter S02). Accounting for the pres- 
ence of point sources was discussed by Hobson et al. (1999) 
for the flat patches and point source detection on the full 
sky maps using Spherical Mexican Hat Wavelets (SMHW) 
was analyzed by Vielva et al. (2003). A joint technique using 
both SMHW and MEM for the flat case was investigated by 
Vielva et al. (2001), and for the spherical case it will be de- 
scribed in a forthcoming paper by Stolyarov et al. (in prepa- 
ration). The method has also been used to construct simu- 
lated all-sky catalogues of the thermal Sunyaev-Zel'dovich 
effect in galaxy clusters (Geisbiisch, Kneissl and Hobson, in 
preparation). 



The separation tests described in previous articles were 
performed making some simplifying assumptions. In partic- 
ular, the receiver noise was assumed uncorrelated and sta- 
tistically homogeneous over the sky, which is a reasonable 
approximation for some scanning strategies. Another simpli- 
fication concerned the foreground spectral behaviour. It was 
assumed that spectral parameters, such as the synchrotron 
spectral index and dust emission properties, were spatially- 
invariant. In the analysis of real data, however, one cannot 
simply ignore these effects, and it is very important to inves- 
tigate their infiuence on the component separation results. 

For real observations, each point in the sky is observed 
a different number of times depending on the scanning strat- 
egy of the instrument. In the case of simple scanning strate- 
gies, such as constant latitude scans, the spin axis stays close 
to the plane of the ecliptic. In this situation, pixels near 
the ecliptic pole are observed several times more often than 
those near the ecliptic plane. This uneven coverage leads to 
marked differences in the noise rms per pixel across the sky. 
We also note that instrinsic gain fluctuations in the instru- 
ment can also contribute to the variable noise rms. 

Component separation using incomplete maps contain- 
ing cuts (e.g. along the Galactic plane) can be considered 
as an extreme case of varying noise rms. One approach is 
to assume that the noise rms for pixels in the cut is for- 
mally infinite (or, equivalently, that they have zero statisti- 
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cal weight). We show below that this method does indeed 
allow the separation method to cope straightforwardly with 
cuts. Moreover, this technique can in principle be applied to 
analyse arbitrarily-shaped regions on the celestial sphere. 

The assumed isotropy of the foreground spectral param- 
eters over the sky is another extreme simplification. For ex- 
ample, the synchrotron spectral index varies in a wide range. 
Giardino et al. (2002) calculated the synchrotron tempera- 
ture spectral index using three low-frequency radio survey 
maps and found it to vary in the range 2.5^ /3408/i420 =S!3.2. 
Variation of the dust colour temperature Tdust is more im- 
portant for the Planck experiment because the HFI chan- 
nels are quite sensitive to thermal dust emission. Schlegel 
et al. (1998) found Tduat to vary in the range 16K to 
20K, about a mean value (rdust) = 18K assuming a single- 
component model. 

The effect of a spatially- varying Galactic dust emissivity 
index /3 on the MEM reconstruction was investigated for case 
of a flat-sky patch by Jones et al. (2000). Several dust sub- 
components with different emissivities, but with the same 
colour temperature, were included in the separation pro- 
cess, which provided a good reconstruction of the compo- 
nents. More recently, Barreiro et al. (2004) used a combined 
real and harmonic space-based MEM technique to perform 
a component separation on real data, in the presence of 
anisotropic noise, cut-sky maps and spectral index uncer- 
tainties. Since this approach requires multiple transitions 
between pixel and spectral domains, however, the computa- 
tion of the necessary spherical harmonic transforms makes 
it is much slower than harmonic-space MEM, and hence it 
was implemented only for low-resolution COBE data. 

In this paper we will demonstrate how to extend 
the full-sky harmonic-space MEM component separation 
method to take into account anisotropic noise and variations 
in spectral parameters, by making use of prior knowledge of 
the uneven sky coverage and the average value of the spec- 
tral parameters. The structure of the paper is as follows. In 
a Section|21we summarise the basics of the MEM component 
separation technique, describe the model of the microwave 
sky used in the simulations, and review the impact of the 
non-isotropic noise on the CMB reconstruction. In a Sec- 
tion El we will introduce an approach for taking account of 
dust colour temperature variations in the component sepa- 
ration process, describe the new microwave sky model with 
variable Tdust and show the results for the reconstruction 
tests made using different approximations. In Section 2] we 
discuss the results and present our conclusions. 



2 ANISOTROPIC NOISE AND INCOMPLETE 
SKY COVERAGE 

Following the notation from previous articles (Hobson 
et al. 1998; S02), a data vector of length Uf, which contains 
the observed temperature (or specific intensity) fluctuations 
at Uf observing frequencies in any given direction x, can be 
defined as 

p Tic 

dy{x) = I Bv{x ■ x)'^^Fyp Sp{x) AQ! + ev{x) (1) 

where F^p is the frequency response matrix, is the beam 
profile for the vth frequency channel, ev{x) is the instru- 



mental noise contribution in the i^th channel, and Sp{x) is 
the signal from pth physical component. An integration is 
performed over the solid angle Q. 

For observations over the whole sky and for a Gaus- 
sian (or at least circularly symmetric) beam profile at each 
frequency we can rewrite the previous equation in matrix 
notation using the spherical harmonic coefficients: 

dem ~ ReStm + etm, (2) 

where dim, a^m and eim are column vectors containing n/, 
ric and n/ complex components respectively. The response 
matrix has dimensions Uf x ric and accommodates the 
beam smearing and the frequency scaling of the components. 

We assume that the anisotropic, uncorrelated pixel 
noise represents a non-stationary random process in the sig- 
nal domain with mean value {e,^{x)) = and an rms that 
varies across the sky. This leads to correlations between dif- 
ferent [l, m) modes in the spherical harmonic domain. As- 
suming that the instrumental noise is uncorrelated between 
different frequency channels, we have 

{eim{ly)el,^,{u')) = Mlm,l'm'{v)S^v' , (3) 

where the form of Mim.tm' {^) is discussed in detail in Ap- 
pendix A. In particular, for the case in which the instrumen- 
tal noise is uncorrelated between pixels (as we consider in 
Section [2. in . we have 

J^em.,i'm' (i^) = f^pix ^ Yi,„{xp)Yi,m' {xp)al{xp), (4) 
p=i 

where a^{xp) — {el{xp)) is the noise variance in the pth 
pixel, f2pix is the pixel area and A^pix is the total number of 
pixels. For each frequency i/, if we define the double indices 
i = £m and j = I'm' , we may consider the quantities I^J as 
the components of the noise covariance matrix M{v) in the 
spherical harmonic domain. It is also shown in Appendix A 
that the elements of the inverse of this matrix are given by 

p=i 

where, in general, the notation [• ■ denotes the ijth ele- 
ment of the corresponding matrix. 

Ideally, one would like to take into account the full noise 
covariance matrix in performing the component separation. 
As discussed in Stolyarov et al. (2002), however, it is not 
computationally feasible to determine the entire vector, a, 
containing the best estimate of the harmonic coefficients 
of the physical components, using all the elements of the 
data vector d simultaneously. Instead, a 'mode-by-mode' ap- 
proach is used, in which any a priori coupling between dif- 
ferent [l, m) modes is neglected. Taking the modes to be in- 
dependent corresponds to assuming that the likelihood and 
prior probability distributions factorise, such that 

Pr(d|a) = J|Pr(d ), (6) 

e,m, 

Pr(a) = ]^Pr(a^™). (7) 

l.m 

This offers an enormous computational advantage, since one 
can maximise the posterior probability 
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Pr(afmjd£m) (X Pr(dft„|afm)Pr(a^„) (8) 

at each mode separately, where Pr(df„i|afm) is the hkehhood 
and Pr(afm) is an entropic prior. 

The factorisation together with the assumption that 
the instrumental noise is Gaussian, implies that the likeli- 
hood function is given by 

Pr(dft„|afm) oc exp [-x^(afm)] , (9) 
where is the standard misfit statistic 

X^i^e-m) = (dfm — Rfafm)^M£m(*^f'" ^ ^e^tm)- (10) 

In this expression is the UfXUf inverse noise covariance 
matrix for the {t, m) mode, which can be different for each 
mode. For the case in which the instrumental noise is uncor- 
related between frequency channels, it is a diagonal matrix 
and we take the i/th diagonal entry to be 

(11) 

where the right-hand side is obtained by setting (! = i 
and m' = m in 0- One should note that we are con- 
sidering inverse noise covariance matrices in the expres- 
sion lllH . This is clearly not equivalent simply to setting 
[Nfm]i^i^ = \M {v')\im,im and then inverting the resulting ma- 
trix (which reduces in this case to taking reciprocals of the 
diagonal elements). 

The reason for using the definition Ijll^ is that it allows 
for the straightforward analysis of cut-sky data. We can con- 
sider the missing areas in cut-sky maps as an extreme case of 
anisotropic noise in which the noise rms of the pixels in the 
cut is formally infinite. As can be seen from (|1J, this leads to 
elements of Miy) which are also formally infinite, and hence 
this causes problems for the analysis. One might try to avoid 
such difficulties by, for example, setting the noise rms in the 
cut to be some arbitrary large value. This can itself cause 
problems of computational accuracy, however, and one also 
has no a priori means of deciding to which large value the 
noise rms in the cut should be set. Moreover, the resulting 
component separation may depend on the value adopted. 
Using @ instead, we see that formally infinite noise rms in 
the cut is easily accommodated. It corresponds simply to 
omitting the pixels in the cut from the summation in @. 
Adopting this approach also means that there is no need 
to apply any smoothing to the edges of the cut, or to per- 
form the analysis on some cut-sky harmonic basis. One is, 
in effect, still performing the component separation over the 
whole sky, but not constraining the solution in the cut re- 
gion in any way. After the component separation has been 
performed one may then apply whichever cut is desired to 
the resulting component maps. 



2.1 Simulated Planck observations 

For our component separation tests, we prepared simu- 
lated Planck observations for 9 frequency bands at 30, 44, 
70 GHz (LFI) and 100, 143, 217, 353, 545 and 857 GHz 
(HFI) with a resolution of 3.4 arcmin (HEALPix parame- 
ter A^side = 1024). The key parameters of the satellite, such 
as beam FWHMs and average noise levels, were taken from 
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Figure 1. The template of the primordial CMB anisotropies used 
in the simulations. The map is plotted in units of /iK 




Figure 2. The noise rms variation across the sky, in Galactic 
coordinates, due to uneven coverage resulting from a sinusoidal- 
precession scanning strategy with 9 = 90°, A6 = 7°, opening 
angle <p = 85° and ricyde = 4. The plot is scaled such that unity 
corresponds to the average noise rms. Note that the noise rms 
decreases near the ecliptic poles. 

the relevant web-page^. In our simulations and separation 
tests, we assumed symmetrical Gaussian beams. 

The models of the astrophysical components used in 
the simulations are the same as those presented in S02. We 
assume the presence of six physical components of emission: 
primordial CMB, thermal and kinetic Sunyaev-Zel'dovich 
(SZ) effects, and Galactic synchrotron, free-free and thermal 
dust emission. Fig. shows the template of the primordial 
CMB anisotropies used in the simulations. 

Anisotropic, uncorrelated instrumental noise is simu- 
lated for the sinusoidal-precession scanning strategy with 
9 = 90°, A9 = 7°, opening angle (f> = 85° and ricydc ~ 4. 
The ecliptic latitude /3 depends on the ecliptic longitude A 
as /3 = A6' sin(ncycioA). A map of the corresponding sky 
coverage for this scanning strategy is shown in Fig. |5| All 
frequency maps are modelled with the same sky coverage for 
simplicity. 

On calculating the 'statistical weight' wern = 
[Af^^]em,em for each {£,m) mode using @, we might expect 
some trends and variations along both £ and m. For this 
particular scanning strategy, however, with full sky cover- 
age the variations are not very large. In Fig. |H]we plot the 
relative amplitude of the statistical weights for the 100 GHz 

^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 
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Figure 3. The statistical weight wim of each {£, m) mode for the 
100 GHz HFI frequency channel assuming the scanning strategy 
illustrated in Fig. 121 In the case of uniform noise ui^,„ = 6.1 X 
lO^'' (MJy/Sr)-2 for this channel. 
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Figure 4. As in Fig. U] but for a map containing a symmetric, 
constant-latitude cut of ±10° along the Galactic plane. 



channel on large scales {£ ^ 200). It is clear from the plot 
that the deviations from the average level do not exceed the 
10 per cent level. 

As mentioned above, a cut-sky map can be considered 
as an extreme case of anisotropic noise. To test the accu- 
racy of the component separation in the presence of a cut, 
we also prepared a set of channel maps with a symmet- 
ric, constant-latitude cut of ±10° applied along the Galac- 
tic plane. In the remainder of the map, we preserved the 
anisotropic noise pattern resulting from the uneven cover- 
age produced by our assumed scanning strategy, but in the 
cut region the noise rms was considered as formally infinite. 
As might be expected, an extreme anisotropy of this sort 
leads to more profound variations in the statistical weight 
wem of the modes. In Fig. |11 we plot wim for the 100 GHz 
channel in this case. It is easy to see from the plot that the 
deviations from the average level are significant. 



2.2 Results for all-sky data 

We first performed a component separation of simulated 
data with the anisotropic noise distribution shown in Fig. |51 
assuming diagonal inverse noise covariance matrices 
with elements given by Ull^ . The reconstructed map ob- 
tained for each physical component was of a similar visual 
quality to those obtained in S02, and so they are not pre- 
sented here. Instead, we plot only the residuals of the recon- 
struction of the primordial CMB component in Fig. |K| We 
see that these residuals are essentially featureless, except for 
a small band in the Galactic plane. We note, however, that 
the residuals do contain a faint imprint of the noise rms 
distribution plotted in Fig. |21 

In Fig. |n] we present the unbiased estimator (see S02) 
of the CMB power spectrum, obtained from the recon- 
structed modes for this component, and compare it with 
the power spectrum of the input primordial CMB realisa- 
tion used in the simulations. We see from the figure that, 
even in the presence of realistic anisotropic noise, the recon- 
structed CMB power spectrum follows the input spectrum 
out to ^ « 2500, recovering the correct positions and heights 
of the first seven acoustic peaks. 

In fact, even without making use of our prior knowledge 
of the noise anisotropy, the effect of the spatially-varying 
noise rms on the reconstructions is negligible. In order to 
compare our new approach with that presented in S02, we 
analysed the same simulated data once more, but assuming 
that the noise was in fact isotropic. The resulting recovered 
maps were visually indistingushable from those obtained us- 
ing the approach outlined above. In Fig|7| we plot the dif- 
ference in the unbiased estimators of the CMB power spec- 
trum obtained by the two methods. We see that the two 
approaches yield almost identical results. This is not too 
surprising, since the variation in the values of the statistical 
weight per {£, m) mode are quite small in the all-sky case, 
as illustrated in Fig. 

2.3 Results for cut-sky data 

In the process of foreground removal from CMB data (as 
opposed to component separation), it is common practice to 
impose some cut along the Galactic plane prior to perform- 
ing the analysis. In S02, it was shown that such an approach 
is unnecessary, since accurate results may be obtained by in- 
stead performing the component separation on all-sky data. 



All-sky component separation 5 




Templ^e PS 

Reconsiruetsd pQwer specirum 

with lo errors 




500 10OO 1500 2000 2500 3000 



Wlultipole number. I 



Figure 8. The reconstruction of the CMB anisotropy in the cut- 
sky case. The map is plotted in units of /iK 



Figure 6. Initial and unbiased reconstructed CMB power spectra 
for the modelled data with non-uniform pixel noise. The Icr errors 
on the reconstructed power spectrum are also shown. 
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Figure 7. The difference in the unbiased estimators of the CMB 
power spectrum obtained assuming uniform noise across the sky, 
and by taking into account noise rms variations. 



and then imposing a Galactic cut in the reconstructed CMB 
map, if required, prior to any further analysis. Nevertheless, 
the question of analysing cut-sky data sets might occur, for 
example, in the event of instrument failure, and it is thus of 
interest to investigate how to accommodate this complica- 
tion. 

To test our approach of modelling cut-sky data as an 
extreme case of anisotropic noise, we performed a compo- 
nent separation using a set of frequency maps with ±10° 
symmetrical, constant-latitude cut along the Galactic plane, 
within which the pixels were set to zero. In Fig. |S|we show 
the reconstruction of the primordial CMB component, and 
the corresponding residuals are shown in Fig.|U] 

The reconstructed physical components were not explic- 
itly restricted to be zero in the cut. Instead, as explained 
above, the reconstruction was simply not constrained by 
the data in the cut region. As a result the reconstruction 
of spherical harmonic modes that lie predominantly in the 



Figure 9. The residuals of the reconstruction in cut-sky case. 
The map is plotted in units of fiK 



cut will be prior driven. Since the maximum of the entropic 
prior occurs for zero signal, such modes will have zero am- 
plitude in the reconstruction. This does not mean, however, 
that the reconstructed temperature map will be precisely 
zero in the cut. Modes that are not restricted to the cut 
region will have amplitudes that are constrained to be non- 
zero by data on the un-cut part of the sky. Such modes may, 
in general, contribute to the reconstructed pixel values in 
the cut region. 

We note from Fig. Inland Fig. El that, by using our ap- 
proach, the quality of the CMB reconstruction outside the 
cut is unaffected by the presence of the cut. In particular, we 
see that the reconstruction residuals are featureless over the 
un-cut part of the sky, and show no identifiable structure 
even directly adjacent to the cut region. In making this ob- 
servation, it should also be remembered that no smoothing 
was applied to the edges of the cut prior to the analysis. 

Finally, in Fig. 1101 we plot the unbiased estimator (see 
S02) of the CMB power spectrum, obtained from the re- 
constructed modes for this component, and compare it with 
the power spectrum of the input primordial CMB realisation 
used in the simulations. We see that, even in the presence 
of the cut, the reconstructed CMB power spectrum follows 
the input spectrum out to £ « 2500. 
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Figure 10. Initial and unbiased reconstructed CMB power spec- 
tra for data containing a Galactic cut. The Icr errors on the re- 
constructed power spectrum are also shown. 



3 SPATIALLY- VARYING FOREGROUND 
SPECTRA 

Another major problem for harmonic-space component sep- 
aration methods, such as MEM and Wiener filtering, is the 
implicit assumption in Q that the emission from each phys- 
ical component can be factorised into a spatial template at 
some reference frequency vq and a frequency dependence, so 
that 



(12) 



This clearly represents the idealised case in which the spec- 
tral parameters of each component do not vary with direc- 
tion on the sky. In terms of the detailed operation of the 
harmonic-space MEM algorithm (and Wiener filtering), this 
assumption is quite central. The MEM approach seeks to 
find the most probable values for the harmonic coefficients 
a of the astrophysical components, given the observed data 
and an entropic regularisation prior on the solution. At each 
step in the optimisation, the predicted data corresponding 
to the current best estimate of a is compared (through the 
Gaussian likelihood function) with the real data. In calcu- 
lating the predicted data at each frequency, one must deter- 
mine the contribution of the physical components at this fre- 
quency by scaling the solution from the reference frequency 
vo using the conversion matrix F^c- This matrix determines 
the frequency behaviour of the components and its coeffi- 
cients are calculated in accordance to some model of spectral 
scaling - power law for the synchrotron, modified blackbody 
emission law for the dust emission, and so on. The scaling 
is taken to be the same for all {£, m) modes, which implies 
that the parameters defining the spectral scaling do not vary 
across the sky. 

The assumption of spatially-invariant spectral parame- 
ters is reasonable for the primordial CMB and the two SZ 
effects. The thermal SZ effect does depend on the electron 
temperature of the clusters Te, which can reach 10-15 keV; 
this dependence is rather weak, but can be used to deter- 
mine Te during component separation, as will be explored 
in a forthcoming paper. For the Galactic components, how- 
ever, our assumption is clearly not valid. For example, as dis- 



cussed in Section the spectral indices of the synchrotron 
and thermal dust emission are known to vary considerably 
across the sky. In the case of Planck observations, variation 
in the synchrotron spectral index will not cause severe prob- 
lems. Even in the lowest frequency channel at i/ — 30 GHz, 
the synchrotron emission is quite weak. A similar conclusion 
may be drawn for free-free emission. Unfortunately, varia- 
tions in the spectral parameters of thermal dust emission 
can give rise to severe difficulties in performing component 
separation for Planck data. 

In the previous simulations presented in S02, the ther- 
mal dust emission was assumed to follow a simple single- 
component grey-body model with spectral emissivity index 
/3 = 2 and mean dust temperature (Tdust) = 18 K. It is 
easy to show that variations in these parameters at 20 — 30 
per cent level can lead to differences of the same factor (in 
specific intensity units) when scaling from the reference fre- 
quency vo = 300 GHz to the highest HFI channel aX u = 857 
GHz. This can cause huge errors in the prediction of dust 
emission which can reach up to 10^ MJy/Sr in some regions 
of the galactic plane. As a result, incorrectly assuming the 
emissivity and dust temperature to be spatially-invariant 
can severely reduce the accuracy with which all the physical 
components are recovered. The dominant cause of these dif- 
ficulties is in fact the uncertainty in the dust temperature, 
which is known to vary in the range 5-25 K. It is there- 
fore necessary to take proper account of the variation in the 
dust temperature to obtain reliable component separation 
results. 

Using the simplest single-component grey-body model 
as an example, the scaling with frequency of the dust emis- 
sion can be defined, in terms of specific intensity, by 



iy°'B{u,T{x)) 
uSB{uo,T{x)) 



I,„ix)F,,{T{x)), (13) 



where Iu{x) is the dust specific intensity at frequency v in 
the direction x, T{x) is dust colour temperature in this di- 
rection, l3 is the emissivity spectral index, which is assumed 
to be uniform over the sky, and B{v,T) is the blackbody 
function. The quantity Fvc{T{x)) is the appropriate element 
of the conversion matrix for the dust component, which de- 
pends on the position on the sky. 

Since M'il takes the form of a product of two functions 
in the spatial domain, we cannot easily pass to the spectral 
domain, because this product will turn into a convolution 
over the whole range of the spherical harmonics. In order 
to work in terms of harmonic coefficients, we instead begin 
by expanding 11311 around the mean dust temperature To to 
obtain 



Iu{x) ^ I,„{x)F,c{To)+I,„{x)AT{x] 



dF,, 



dT 



+ ■ 



T=Tq 



(14) 



where AT{x) is the deviation from the mean temperature in 
a given direction. We have not included further terms in the 
expansion for the sake of brevity, but these may be written 
down trivially. Given our Taylor series expansion, it is now 
possible to use our standard formalism to reconstruct sev- 
eral separate (but highly correlated) fields - intensity at the 
reference frequency I^q and intensity-weighted temperature 
fields IiygAT, /yo(Ar)^, and so on, to obtain an increasingly 
accurate approximation to Iv{x). 



All-sky component separation 7 



A similar approach may also be used for accommodating 
spatial variation in any other foreground spectral parameter, 
such as the dust emissivity or synchrotron spectral index. 
One may also treat the electron temperature of SZ clusters in 
the same way. It should be noted, however, that the number 
of terms used in the Taylor expansion cannot simply increase 
without limit. Typically, one is constrained such that the 
total number of fields to be reconstructed does not exceed 
the number of frequencies at which observations are made. 

3.1 Simulated Planck observations 

The simulated Planck observations considered here are 
similar to those described in Section 12.11 We assume six 
physical components and simulate nine frequency maps, 
smoothed with appropriate FWHMs and add anisotropic 
pixel noise. The only difference here is that the thermal dust 
emission has a spatially-varying temperature distribution. 

In simulating the thermal dust emission, it is more ac- 
curate to use a two-component best-fit thermal dust model, 
with the dominant component having a mean temperature 
of (Ti) — 16.2 K and a small contribution of cold dust with 
{T2) = 9.4 K (Finkbeiner et al. 1999). The two compo- 
nents are fully correlated and their temperature is related by 
T2 = 0.352 Ti^'^*, so we have only one effective temperature 
parameter. A template for the temperature Ti{x) across the 
sky was constructed from the colour temperature of Schlegel 
et al. (1998). The latter was constructed assuming a single- 
component dust model with the uniform emissivity index 
of /3 = 2 and has the mean temperature (Tdust) = 18K. 
This map was simply scaled to (Tdust) = 16. 27^ to obtain 
the T\[x) map for the two-component model used in our 
simulations (see Fig. IIH . 

For a multi-component dust model, the dust emission, 
in specific intensity units, scales with frequency as 



= /.„(4)F,.,(r(4)), 



(15) 



where k = 1,2,... labels the dust component, fk is a nor- 
malisation factor for the fcth component and fk = 1- 
The terms qkiy/vo)"'' represent the relative emission effi- 
ciency for each component (see Finkbeiner et al (1999) for 
details). For the two-component best-fit model used in our 
simulations, /i =0.0363, q:i=1.67, 0:2=2. 70 and gi/g2 = 13.0 
. We can expand the multi-component dust model around 
the mean temperature an analogous manner to that given 
in 11411 for the single component case. 



3.2 Component separation results 

Several separation tests were performed to investigate the 
approach outlined above. In our first test, we recreated the 
method used in S02, in which the dust temperature was 
assumed to be spatially-invariant (the zeroth-order approxi- 
mation), and so only intensity fields ll!^ (A; = 1, 2, . . . , 6) for 
the six physical components were reconstructed at the refer- 
ence frequency vq. In the second test, we attempted to recon- 
struct 7 components, which included 11^ for = 1, 2, . . . , 6, 
but also the intensity-weighted dust colour temperature 
field I^gAT (the first order of the expansion). Finally, in a 




Figure 11. The map of the temperature variations across the 
sky in K (Schlegel et al. 1998). The original map from Schlegel 
et al. 1998, with {T^^st) = ^8K, was scaled to (Tdust) = 16.2^" 
and the mean value has been subtracted. 







Figure 12. The reconstruction of the colour dust temperature 
variations ATdust- The map is plotted in units of K 



third test, 8 components were reconstructed, namely I^p' for 



fc = 1, 2, . . . , 6, together with hgAT and 7^(,(AT)^, thereby 
accounting for dust temperature variations up to second or- 
der. 



3.2.1 Quality of the CMB reconstruction 

The reconstruction of the primordial CMB component and 
the corresponding residuals are shown in Fig. llSl for all three 
separation tests. It is clear from the figure that the separa- 
tion assuming Tdust to be constant over the sky has a very 
strong residual signal along the plane of Galaxy on the CMB 
map, where the dust component is strong. Moreover, we see 
that the reconstruction outside the Galactic plane is also 
badly contaminated; this is discussed further in Section 0| 
We thus conclude that the practice adopted in previous com- 
ponent separation analyses of assuming a spatially-invariant 
dust temperature leads to CMB reconstructions of unac- 
ceptably poor quality. Although we do not plot them here, 
the reconstructions of the other physical components are 
found to contain similar artefacts both inside and outside 
the Galactic plane. 

Fortunately, the reconstruction errors are significantly 
reduced even when the dust temperature variations are ac- 
counted for only up to first order in the expansion 11141 . We 
see that outside Galactic regions of high dust emission, the 
CMB residuals in this case are essentially featureless. Never- 
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Figure 13. Reconstructions of the primordial CMB component (left-hand column) and the corresponding residuals (right-hand column) 
for the three separation tests: without any accounting for Tjust variations across the sky (top panel); accounting for variations up to first 
order (middle panel); and accounting for variations up to second order (bottom panel). All maps are plotted in units of ^K. 



theless, the Galactic plane still contains significant dust con- 
tamination. Tiie quality of tiie CMB reconstruction can be 
improved still furtiier by accounting for Tdust variations up 
to second order. In this case, the remaining contamination 
of the CMB residuals in the Galactic plane is significantly 
reduced, to a level comparable to that obtained in S02 in 
the analysis of simulated data in which Tdust was taken not 
to vary. 

A more quantitative comparison of the quality of the 
CMB reconstructions is given in Fig. ll4l (left panel), in which 
the unbiassed estimates (see S02) of the angular power spec- 
tra of the three reconstructions are plotted in relation to the 
input power spectrum. In each case, the power spectrum has 
been calculated using the full sky reconstruction. As antic- 
ipated from Fig. ^| the recovered power spectrum in the 
zeroth-order approximation is very poor due to the signif- 
icant dust contamination. The recovered power spectrum 
in the first-order approximation is significantly improved. 



It follows the input spectrum out to i 1700, recovering 
the positions and heights of the first 5 acoustic peaks. Be- 
yond this point, the presence of the 6th and 7th acoustic 
peaks is inferred, but the recovered spectrum begins to un- 
derestimate the true power level. In the second-order ap- 
proximation, the power spectrum recovery is improved still 
further, and the positions and heights of the first 7 acous- 
tic peaks are accurately recovered before the true power 
level is again underestimated. The differences between the 
first and second-order recovered power spectra are further 
illustrated in the right panel of Fig. 1141 in which the ratio 
{Cr - cf'"P'"*'=)/C7f""P'"*° is plotted for these two cases. 
The improved accuracy and extended range of validity of 
the second-order approximation is easily seen. 



All-sky component separation 9 




3.2.2 Reconstruction of the dust temperature distribution 

For the second and third separation tests, the actual dust 
colour temperature variations can be calculated as 

Tdust(i) = (Tdust) + {h,AT)T'/ltr\^)^ (16) 

where 1^"''* is the reconstruction of dust intensity at the 
reference frequency vo assuming the average colour tem- 
perature (Tdust), and (JynAr)'^"^* is the reconstruction of 
intensity-weighted temperature variation field. The result 
is shown in Fig. 1121 Comparing with Fig. 1111 we see that 
the reconstuction has faithfully recovered the main features 
of the input map. It is clear, however, that ATduat varia- 
tions are better restored around the Galactic plane. This 
is not surprising, since we have actually reconstructed the 
intensity- weighted Tdust field, which can be recoverd more 
accurately in areas with higher dust intensity, which are ob- 
viously located predominantly in the Galactic plane. 



4 DISCUSSION AND CONCLUSIONS 

In this paper, we have demonstrated an approach that al- 
lows one to account for the spatial variations of the noise 
properties and spectral characteristics of foregrounds in 
the harmonic-space maximum-entropy component separa- 
tion technique. 

In Section 12.21 we show that the impact of a realistic 
level of the anisotropic noise on the quality of the foreground 
separation is quite small, at least for the simple scanning 
strategy assumed in generating our simulated observations. 
This can be explained by the fact that the rms noise level 
differs from the average level only in relatively small areas 
around ecliptic poles. This leads to variations in the noise 
level at each harmonic mode that differ from the average 
value by only a few per cent. We also illustrate in Section r2.2l 
that the it is possible to perform harmonic-space component 
separation on cut-sky maps by treating the cut as an ex- 
treme example of anisotropic noise. This appraoch has the 
advantage of not requiring one to smooth the edges of the 
cut with some apodising function prior to the analysis. 



In Section |3] we show that the variation of spectral pa- 
rameters of foregrounds may be taken into account by a 
method of succesive approximations based on a series ex- 
pansion of the corresponding intensity field around the mean 
value of the parameter. In particular, we investigate the ef- 
fect of dust colour temperature variations on the quality 
of the component separation, focussing on the reconstruc- 
tion of the CMB. We show that realistic dust temperature 
variations lead to severe contaminaton of the CMB recon- 
struction if, in the separation process, the dust temperature 
is assumed not to vary. This contamination is concentrated 
in the Galactic plane, but significant artefacts exist at high 
Galactic latitudes. The poor quality of the reconstruction 
outside the Galactic plane is a result of performing the re- 
construction mode-by mode in harmonic space. The inaccu- 
rate model of the dust emission leads to errors in the de- 
termined amplitudes of a wide range of spherical harmonics 
in the CMB reconstruction. Many of these modes do not lie 
predominatly in the Galactic plane region, but contribute to 
the reconstruction over the whole sky. 

If one is content simply with removing foregrounds from 
the CMB, rather than performing a component separation, 
then one could apply a Galactic cut prior to the analysis. 
In Fig. 1151 we show the CMB reconstruction and residuals 
respectively obtained by applying a Galactic cut of ±25° 
to each simulated frequency map used in Section 13.11 and 
assuming a constant dust temperature across the sky. The 
noise rms in the cut region was assumed to be formally infi- 
nite in the manner discussed in Section |5| We see from the 
figures that the quality of the reconstruction is significantly 
improved as compared with the case in which no Galactic 
cut was applied, which was shown in Fig. 1131 (top row). In 
particular, we note that the residuals now contain no obvi- 
ous artefacts outside of the cut region. Thus, even assuming 
an inaccurate dust model, one can still recover a reasonable 
reconstruction of the CMB outside of the Galactic plane. 

It is clearly not possible, however, to perform an accept- 
able all-sky component separation for the Planck experi- 
ment by assuming constant dust spectral parameters, even 
using all 9 frequency channels. Nevertheless, taking account 
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Figure 15. The reconstruction (left) and residuals (right) of the CMB for the case plotted in Fig. 1131 (upper row), but with a Galactic 
cut of ±25°. The maps are plotted in units of fiK. 



of dust temperature variations up to first order in the series 
expansion significantly improves the CMB reconstruction to 
an acceptable level. This reconstruction quality is still fur- 
ther improved by including second-order terms and is then 
comparable to that obtained for the ideal case presented in 
Stolyarov et al. (2002), where the simulated observations as- 
sumed no dust temperature variation across the sky. More- 
over, an accurate reconstruction of the dust temperature 
variation is obtained over the whole sky. 

Finally we note that the approach for dealing with 
spatially-varying spectral parameters described in this paper 
can also be applied to other foregrounds to yield, for exam- 
ple, maps of the synchrotron spectral index. The method 
can also be used to reconstruct the electron temperature in 
clusters from their thermal SZ effect, as will be discussed in 
a forthcoming paper. 
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APPENDIX A: THE NOISE COVARIANCE 
MATRIX IN THE SPHERICAL HARMONIC 
DOMAIN 

We consider a single frequency map for which rip is the in- 
strumental noise in the pth pixel. We denote the spherical 
harmonic coefficients of the noise by etm- 

The elements of the noise covariance matrix in the har- 
monic domain are given by 

p p' 
V p' 

= fipix^i^fm(p)lf*m'(p)crp (Al) 



All-sky component separation 



where yim(p) denotes the value of the corresponding spheri- 
cal harmonic at the pth pixel centre, fipix — 4ir/Npi^ and, in 
the last line, we have assumed that the noise is uncorrelated 
between pixels, so that {upUpi) = apSppi . We note that in 
the special case in which the noise is statistically isotropic. 



a for all pixels, and we obtain 



Afe, 



S^pix ^ Ye„,{p)Yt1^, (p) 



(A2) 



It is convenient to define the double indices i = £m 
and j = £'m' , and regard ijAljl as the elements of the noise 
covariance matrix Af. One may then write 



A/' = npixY'''NY, 



(A3) 



where N is the noise covariance matrix in the pixel domain, 
with elements Nppi = ijipUpi), and we have defined the 
transformation matrix Y with elements Ypi = Yiijp). 

Using the fact that transformation matrix has the useful 
property 



YTy 



1 

^pix 



I = 



(A4) 



it is straightforward to show that AfAf^^ — I = Af^^Af, 
where the inverse noise covariance matrix in the spherical 
harmonic domain is given by 

= ytN^W. (A5) 

In the special case that the instrumental noise is uncorre- 
lated between pixels, N^^ is diagonal and we obtain 



Ytm (P) Vf'm' (P) • 



(A6) 



This paper has been typeset from a Tff)i/ IW^^X file prepared 
by the author. 



